Control estadístico de procesos
Introducción a R
Jordi Cuadros, Lucinio González
Diciembre de 2017

Elementos básicos de R

R

R es un entorno y un lenguaje de programación estadística.

Es una de las herramientas más usadas en el análisis de datos. Es de código libre, su uso es gratuito y dispone de una comunidad muy activa que contribuye al desarrollo de herramientas y funciones adicionales (más de 10.000 paquetes) y a la resolución de preguntas a través de la red.

https://www.r-project.org/about.html https://www.rdocumentation.org/ https://stackoverflow.com/questions/tagged/r

Instalación de R y RStudio

R como herramienta de cálculo

R puede usarse desde la consola introduciendo directamente las instrucciones correspondientes.

2 + 3 
## [1] 5
8 ^ 10
## [1] 1073741824
log(100) # Logaritmo neperiano
## [1] 4.60517

# indica que lo que sigue es un comentario

Uso de scripts

Un script en R es un conjunto de instrucciones escritas en un archivo de texto, de forma que estás puedan almacenarse y ejecutarse con posterioridad. Se suelen almacenar con la extensión .R.

Un editor muy recomendado para la creación, ejecución y depuración de estos scripts es RStudio.

Acceso a la ayuda

  • Buscar ayuda sobre una función
? "mean"
  • Buscar un texto en la ayuda
?? "anova"

En RStudio, se puede seleccionar y pulsar F1 para buscar en la ayuda

Datos básicos

  • numeric: número decimal de 15 dígitos significativos
  • integer: entero de 32 bits (hasta 2·109 aprox.)
  • logical: TRUE o FALSE
  • character: cadena de caracteres de longitud indeterminada

Datos básicos – numeric

a <- 2
a
## [1] 2
class(a)
## [1] "numeric"
b <- 13.6788956789
b
## [1] 13.6789
class(b)
## [1] "numeric"
print(b, digits = 10)
## [1] 13.67889568

Operadores más comunes para datos de tipo numeric

Suma +
Resta -
Producto *
División /
Potencia ^ o **
Módulo (residuo) %%
División entera %/%
Comparaciones == > < >= <= != 
a+b
## [1] 15.6789
a-b
## [1] -11.6789
a*b
## [1] 27.35779
a^b
## [1] 13114.69
a/b
## [1] 0.1462106

Datos básicos – integer

Permite almacenar números enteros. Es el tipo de datos usado para contadores e índices.

n <- as.integer(340000)
class(n)
## [1] "integer"
n <- 2L
class(n)
## [1] "integer"

Datos básicos – character

Permite almacenar cadenas de texto de longitud indeterminada.

Los literales se indican entre comillas dobles (o simples).

a <- "aaa"
b <- "bbb"

paste(a, b, "hola", sep = ", ") 
## [1] "aaa, bbb, hola"

Operadores más comunes para datos de tipo character

Comparaciones == > < >= <= !=

Datos básicos – logical

Almacena valores que son verdadero (TRUE) o falso (FALSE).

Es el resultado de una comparación.

3 == 2
## [1] FALSE
b <- 3 != 2 
b
## [1] TRUE

Operadores más comunes para datos de tipo logical

Intersección lógica &
Unión lógica |
Negación lógica !
Comparaciones == !=

Por ejemplo…

a <- TRUE
b <- F
a & b # Operador AND
## [1] FALSE
a | b # Operador OR
## [1] TRUE
!b  # Operador NOT
## [1] TRUE

Conversión a datos básicos

Las funciones de conversión de tipos en R se indican como as. seguido del tipo de dato de destino.

as.character(2)
## [1] "2"
as.numeric(TRUE)
## [1] 1

Para saber si una información o variable es de un tipo determinado puede usarse is. seguido del tipo de dato a comprobar.

Constantes

R incorpora algunas constantes preestablecidas como

pi
## [1] 3.141593
Inf
NaN
NA
NULL

También existen funciones para evaluar si un dato corresponde a una de estas constantes.

is.na(3)
## [1] FALSE
is.null(NULL)
## [1] TRUE
is.infinite(-Inf)
## [1] TRUE

Datos compuestos

  • Principales
    • Vector
    • Factor
    • Data frame (Tabla de datos)
  • Otros tipos comunes
    • Factor ordenado
    • Lista
    • Matriz

Datos compuestos – vector

a <- c(2, 3, 4)
str(a)
##  num [1:3] 2 3 4
a[2]
## [1] 3

Las operaciones se aplican a vectores y devuelven vectores

a + a 
## [1] 4 6 8
a > 2.5
## [1] FALSE  TRUE  TRUE

Atención al reciclaje de datos…

a <- c(2, 3, 4)
b <- c(10, 20)
a * b
## Warning in a * b: longer object length is not a multiple of shorter object
## length
## [1] 20 60 40

Las funciones en R suelen ser vectoriales

  • algunas devuelven vectores
a <- c(2, 3, 4)
abs(sin(a))
## [1] 0.9092974 0.1411200 0.7568025
exp(a)
## [1]  7.389056 20.085537 54.598150
  • otras devuelven valores agregados
length(a)
## [1] 3
sum(a)
## [1] 9
mean(a)
## [1] 3

También sd, max, min…

Puede determinarse si un valor está presente en un vector con el operador %in%.

a <- c(2, 3, 4)
3 %in% a
## [1] TRUE
c(2,5,3,4) %in% a
## [1]  TRUE FALSE  TRUE  TRUE

Existen estructuras especificas para crear vectores secuencia.

1:10
##  [1]  1  2  3  4  5  6  7  8  9 10
seq(1, 10, by = 2)
## [1] 1 3 5 7 9
seq(1, 3, length.out = 5) 
## [1] 1.0 1.5 2.0 2.5 3.0

Los vectores se ordenan usando las funciones sort y order.

a <- c(8, 2, 5, 3)
sort(a)
## [1] 2 3 5 8
a[order(a)]
## [1] 2 3 5 8
a[order(a,decreasing = T)]
## [1] 8 5 3 2

Datos compuestos – factor

Un factor es un vector de cadenas indexado. Normalmente se crea a partir del vector de cadenas de texto.

a <- c("hola", "adeu","hola", "adeu", "adeu", "bye")
b <- as.factor(a)
as.character(b)
## [1] "hola" "adeu" "hola" "adeu" "adeu" "bye"
as.numeric(b)
## [1] 3 1 3 1 1 2
str(b)
##  Factor w/ 3 levels "adeu","bye","hola": 3 1 3 1 1 2

La función factor permite la codificación o recodificación manual de un vector.

a <- factor(c(3, 1, 3, 1, 1, 2), labels = c("adeu", "bye", "hola"))
a
## [1] hola adeu hola adeu adeu bye 
## Levels: adeu bye hola
levels(a)
## [1] "adeu" "bye"  "hola"

Datos compuestos – data frame

dfA <- data.frame(int = 1:10,
                  let = sample(letters, 10, replace = TRUE), 
                  ran = rnorm(10))
dfA
##    int let        ran
## 1    1   w  0.3434244
## 2    2   q  0.2019293
## 3    3   i  1.6396732
## 4    4   f  1.8508826
## 5    5   h  0.3068050
## 6    6   n  0.2221603
## 7    7   o  1.8145489
## 8    8   f -0.5816948
## 9    9   e -1.0427837
## 10  10   o -0.7087821
dim(dfA) # Dimensión
## [1] 10  3
nrow(dfA) # Número de columnas
## [1] 10
ncol(dfA) # Número de filas
## [1] 3
str(dfA)
## 'data.frame':    10 obs. of  3 variables:
##  $ int: int  1 2 3 4 5 6 7 8 9 10
##  $ let: Factor w/ 8 levels "e","f","h","i",..: 8 7 4 2 3 5 6 2 1 6
##  $ ran: num  0.343 0.202 1.64 1.851 0.307 ...
head(dfA, 3) # Primeros datos, 6 por defecto
##   int let       ran
## 1   1   w 0.3434244
## 2   2   q 0.2019293
## 3   3   i 1.6396732
tail(dfA, 2) # Últimos valores
##    int let        ran
## 9    9   e -1.0427837
## 10  10   o -0.7087821

Para acceder a los datos almacenados en el data frame se usan índices. Las variables de data frame también pueden extraerse usando su nombre.

dfA[2,3]
## [1] 0.2019293
dfA[,1]
##  [1]  1  2  3  4  5  6  7  8  9 10
dfA$let
##  [1] w q i f h n o f e o
## Levels: e f h i n o q w

Veremos más formas de acceder y usar los datos de una tabla de datos más adelante.

Más datos compuestos – lista

a <- list(2, "2", FALSE)
b <- list(3, "hola", c(2, 3, 4))
a
## [[1]]
## [1] 2
## 
## [[2]]
## [1] "2"
## 
## [[3]]
## [1] FALSE
b
## [[1]]
## [1] 3
## 
## [[2]]
## [1] "hola"
## 
## [[3]]
## [1] 2 3 4
length(a)
## [1] 3
a[[3]]
## [1] FALSE
b[[3]][1]
## [1] 2
  

str(b)
## List of 3
##  $ : num 3
##  $ : chr "hola"
##  $ : num [1:3] 2 3 4

Más datos compuestos – factor ordenado

notes <- c("Aprovat", "Insuficient", "Notable", "Insuficient",
           "Notable", "Excel·lent", "Aprovat")
notes <- factor(notes,
        levels = c("Insuficient", "Aprovat",
                   "Notable", "Excel·lent"),
        ordered = TRUE)
str(notes)
##  Ord.factor w/ 4 levels "Insuficient"<..: 2 1 3 1 3 4 2
levels(notes)
## [1] "Insuficient" "Aprovat"     "Notable"     "Excel·lent"

Más datos compuestos – matriz

a <- matrix(c(2, 4, -3, 5), ncol = 2)
a
##      [,1] [,2]
## [1,]    2   -3
## [2,]    4    5
a[2,2]
## [1] 5
a * a # Producto posición por posición
##      [,1] [,2]
## [1,]    4    9
## [2,]   16   25
a %*% a # Producto matricial
##      [,1] [,2]
## [1,]   -8  -21
## [2,]   28   13
t(a) # Transposición
##      [,1] [,2]
## [1,]    2    4
## [2,]   -3    5

Instalación y carga de paquetes en R

R tiene muchos paquetes para resolver problemas específicos. Para usar un paquete adicional este debe instalarse y cargarse en memoria.

installed.packages()[,1] # Lista los paquetes instalados
(.packages()) # Lista los paquetes en memoria

Dos recursos importantes para buscar y identificar paquetes relevantes son https://www.rdocumentation.org/ https://cran.rstudio.com/web/views/

Para instalar un paquete, por ejemplo “nycflights13”

install.packages("nycflights13")

Para cargar un paquete en memoria, se usa

library("nycflights13")

En RStudio la gestión de paquetes también puede hacerse desde la interfaz del programa

En un script y para garantizar la disponibilidad de un paquete

if(!require("nycflights13")) {
  install.packages("nycflights13")
  library("nycflights13")
}

Para acceder a la documentación de un paquete

help(package="nycflights13")

Algunos paquetes tienen información adicional a la que se puede acceder con las funciones vignette() y demo().

Gráficos en R

¿Para qué usamos los gráficos?

La principales funciones de los gráficos en el análisis de datos son

  • explorar los datos y facilitar su comprensión,
  • permitir el descubrimiento de patrones no evidentes, y
  • comunicar efectivamente los resultados de los análisis.

Gráficos en R

Veremos dos paradigmas distintos para la creación de gráficos en R.

  • Gráficos creados con el paquete base.
  • Gráficos creados en base a la gramática de gráficos (GoG), usando el paquete ggplot2

Gráficos en base R

Los gráficos en base R se construyen a partir de datos en vectores y usando funciones específicas en función del gráfico a realizar.

Para ejemplificar algunas de estas funciones, usaremos el conjunto de datos mtcars, que forma de los paquetes básicos de R.

str(mtcars)
## 'data.frame':    32 obs. of  11 variables:
##  $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
##  $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
##  $ disp: num  160 160 108 258 360 ...
##  $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
##  $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
##  $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
##  $ qsec: num  16.5 17 18.6 19.4 17 ...
##  $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
##  $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
##  $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
##  $ carb: num  4 4 1 1 2 1 4 2 2 4 ...
head(mtcars, 10)
##                    mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
## Duster 360        14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
## Merc 240D         24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
## Merc 230          22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
## Merc 280          19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4

Veremos como crear

  • gráficos de dispersión,
  • histogramas,
  • diagramas de barras, y
  • diagramas de caja (boxplot)

Gráficos en base R – gráfico de dispersión

plot(mtcars$disp,mtcars$hp)

Gráficos en base R – histograma

hist(mtcars$disp)

Añadiendo una curva de densidad…

hist(mtcars$disp,breaks=20,freq=FALSE)
lines(density(mtcars$disp))

Gráficos en base R – diagrama de barras

barplot(table(as.factor(mtcars$cyl)))

Gráficos en base R – diagrama de caja

boxplot(mtcars$disp)
points(mean(mtcars$disp))

Gramática de gráficos (GoG)

La gramática de gráficos es una aproximación teórica al estudio de los componentes de un gráfico. De acuerdo con este análisis, un gráfico se puede construir mediante la especificación de un conjunto de capas y componentes que definen los datos, la asociación de los mismos a aspectos del gráfico, la especificación de la relación entre los valores de las variables de los datos con las del gráfico, la estructura geométrica del gráfico…

Wilkinson, L. (2006). The grammar of graphics. Springer Science & Business Media.

ggplot2

ggplot2 es un paquete de R que implementa de la gramática de gráficos.

Wickham, H. (2010). A layered grammar of graphics. Journal of Computational and Graphical Statistics, 19(1), 3-28.

Referencias:

ggplot2 forma parte del paquete tidyverse (aunque también puede instalarse y cargarse autónomamente).

if(!require("tidyverse")) {
  install.packages("tidyverse")
  library("tidyverse")
}

ggplot2 – capas y elementos del gráfico

En ggplot2 cada elemento gráfico que representa un conjunto de datos constituye una capa. Una o más capas constituyen un gráfico.

Cada capa queda definida mediante la especificación de sus elementos. Los principales son

  • datos,
  • mapeado estético, y
  • geometrías

En ggplot2, los gráficos constituyen un objeto de R y se construyen de forma aditiva.

Por ejemplo,

grafico <- ggplot(data = anscombe,
        mapping = aes(x = x1, y = y1))  # Datos y mapeado estético 
grafico <- grafico + geom_point()       # Geometría

grafico

ggplot2 – datos

En ggplot2, el elemento data (datos) se introduce como primer argumento de la función ggplot. Debe corresponder a una tabla de datos o un tipo de datos convertible a tabla de datos.

grafico <- ggplot(data = anscombe,

ggplot2 – mapeado estético

El mapping (mapeado estético) corresponde al establecimiento de relaciones entre variables de los datos y variables del gráfico. Es el segundo argumento de la función ggploty debe crearse con al función de apoyo aes.

        mapping = aes(x = x1, y = y1))

Para variables cuantitativas, los mapeados más comunes corresponden a

  • posiciones: x, y…
  • tamaño: size
  • color: color, fill

Para variable cualitativas, los mapeados más frecuentes son

  • posiciones: x, y…
  • color: color, fill
  • forma: shape

ggplot2 – geometrías

Las geometrías (geom_) indican la forma que debe tener el gráfico, es decir, cómo se articulan las variables del gráfico. Se añaden al gráfico sumándose al objeto creado por ggplot.

grafico <- grafico + geom_point()

Son geometrías de uso común

  • geom_point
  • geom_line, geom_vline, geom_hline
  • geom_bar
  • geom_histogram
  • geom_boxplot

Un resumen de las geometrías y su relación con las variables del gráfico que reconoce cada una de ellas figura en https://github.com/rstudio/cheatsheets/raw/master/data-visualization-2.1.pdf.

Gráficos en ggplot2

Veremos cómo crear en ggplot2 los gráficos más habituales, añadiendo algunas consideraciones para aquellos casos donde los gráficos realizados con base tienen prestaciones insuficientes.

  • gráficos de dispersión,
  • histogramas,
  • diagramas de barras, y
  • diagramas de caja (boxplot).

De forma general, los gráficos en ggplot2 se construyen a partir de tablas de datos (data frames), de los cuales se seleccionan las variables a representar.

Usaremos 1000 datos del conjunto de datos diamonds para crear los distintos ejemplos.

diaM <- diamonds[sample(1:nrow(diamonds),1000),]
str(diaM)
## Classes 'tbl_df', 'tbl' and 'data.frame':    1000 obs. of  10 variables:
##  $ carat  : num  0.41 1.16 0.3 0.41 0.24 1.11 0.51 1.31 0.41 1.16 ...
##  $ cut    : Ord.factor w/ 5 levels "Fair"<"Good"<..: 2 4 4 5 3 5 5 5 5 5 ...
##  $ color  : Ord.factor w/ 7 levels "D"<"E"<"F"<"G"<..: 2 4 1 5 2 4 4 4 2 6 ...
##  $ clarity: Ord.factor w/ 8 levels "I1"<"SI2"<"SI1"<..: 3 4 6 8 6 8 6 4 5 2 ...
##  $ depth  : num  63.5 59.9 62.4 62.9 61.6 61.2 61.8 61.9 62.2 62.4 ...
##  $ table  : num  57 60 58 54 59 54 56 56 55 55 ...
##  $ price  : int  755 7203 965 1243 552 10053 1875 9085 1030 3800 ...
##  $ x      : num  4.69 6.88 4.3 4.74 3.98 6.71 5.1 7.08 4.79 6.76 ...
##  $ y      : num  4.73 6.85 4.26 4.71 4.01 6.73 5.16 7.05 4.76 6.69 ...
##  $ z      : num  2.99 4.11 2.67 2.97 2.46 4.11 3.17 4.37 2.97 4.2 ...

Gráficos en ggplot2 – gráfico de dispersión

ggplot(diaM, aes(x=carat,y=price)) + geom_point()

Añadiendo una tercera variable (cut) y modificando algunos aspectos de formato…

ggplot(diaM, aes(x=carat,y=price,color=cut)) + 
  geom_point(alpha=.8,shape=21,size=3)

Añadiendo líneas de tendencia…

ggplot(diaM, aes(x=carat,y=price,color=cut)) + 
  geom_point(alpha=.8,shape=21,size=3) +
  geom_smooth(method="lm",se=FALSE)

Gráficos en ggplot2 – histograma

ggplot(diaM, aes(x=price)) + geom_histogram(binwidth=1000)

Y en función del corte…

ggplot(diaM, aes(x=price,fill=cut)) +
  geom_histogram(position='dodge',binwidth=1000)

En frecuencias relativas (por grupo)…

ggplot(diaM, aes(x=price,y=..density..,fill=cut)) +
  geom_histogram(position='dodge',binwidth=1000)

Quizás funcione mejor un gráfico de densidades…

ggplot(diaM, aes(x=price,fill=cut)) +
  geom_density(alpha=.3)

Gráficos en ggplot2 – diagrama de barras

ggplot(diaM, aes(x=clarity)) + geom_bar()

En función de la claridad…

ggplot(diaM, aes(x=clarity, fill=cut)) + geom_bar()

Para comparar entre frecuencias absolutas, funcionan mejor las barras separadas.

ggplot(diaM, aes(x=clarity, fill=cut)) + geom_bar(position="dodge")

Para comparar entre frecuencias relativas acumuladas, son mejores las barras apiladas en frecuencia relativa (para cada clase).

ggplot(diaM, aes(x=clarity, fill=cut)) + geom_bar(position="fill")

Los diagramas de barras también pueden crearse a partir de tablas de datos agrupados. En este caso, debe indicarse qué variable es la y e incluir stat="identity" en el geom_bar.

Gráficos en ggplot2 – diagrama de caja

ggplot(diaM, aes(x=1, y=price)) + geom_boxplot()

NOTA: Para el geom_boxplot se requieren dos variables. Si no hay variable independiente se puede incluir x=1.

Y en función del corte…

ggplot(diaM, aes(x=cut, y=price)) + geom_boxplot()

El gráfico se puede mejorar mostrando todos los puntos, con una posición aleatorizada.

ggplot(diaM, aes(x=cut, y=price)) + 
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(shape = 21, alpha=.5,height=0,width=.2)

O incluyendo un violin plot y un punto para la media…

medias <- diaM %>% group_by(cut) %>%
  summarise(price=mean(price))

ggplot(diaM, aes(x=cut, y=price)) + 
  geom_violin() +
  geom_boxplot(outlier.shape = NA, width = 0.1) +
  geom_point(data=medias,shape=3)

ggplot2 – elementos adicionales

Además de los elementos ya presentados (datos, mapeado estético y geometrías), otros elementos de ggplot2 permiten controlar aspectos adicionales del gráfico, por ejemplo

  • scale_...: controlan los aspectos relativos a la presentación de las variables del gráfico,
  • coord_...: establecen el sistema de coordenadas usado en la geometría,
  • labs: establece los títulos del gráfico,
  • theme, theme_...: controlan la part del gráfico que no corresponde a datos (non-data ink), y
  • facet_: permite la creació de secuencia de gráficos en función de una o dos variables

Un ejemplo para terminar…

ggplot(diaM, aes(x=carat, y = price, shape = cut, col = clarity)) +
  geom_point(alpha=.6) +
  scale_x_continuous(breaks=1:3) +
  scale_y_continuous(trans="log10") +
  scale_color_brewer(palette="Spectral")+
  facet_grid(cut ~ clarity) + 
  theme_bw() + 
  theme(legend.position = "none",text = element_text(size=10))

Estadística básica con R

Funciones de la estadística

La estadística, como parte de la matemática que se ocupa de la recolección, análisis e interpretación de datos, tiene dos funciones principales

  1. la descripción de conjuntos de datos (estadística descriptiva), y
  2. la extracción de conclusiones a partir de estos datos (inferencia)

En este repaso de la estadística básica, usaremos (de nuevo) el conjunto de datos mtcars.

str(mtcars)
## 'data.frame':    32 obs. of  11 variables:
##  $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
##  $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
##  $ disp: num  160 160 108 258 360 ...
##  $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
##  $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
##  $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
##  $ qsec: num  16.5 17 18.6 19.4 17 ...
##  $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
##  $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
##  $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
##  $ carb: num  4 4 1 1 2 1 4 2 2 4 ...
head(mtcars, 10)
##                    mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
## Duster 360        14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
## Merc 240D         24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
## Merc 230          22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
## Merc 280          19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4

Descriptiva de una variable

A menudo, el primer interés que tendremos es el de describir un conjunto de datos correspondientes a una sola variable. Esta es la situación con la que empezaremos.

¿Cómo se describe un conjunto de datos univariante?

Normalmente se consideran los siguientes aspectos

  • frecuencias y cuantiles,
  • tendencia central (también llamado localización o posición),
  • dispersión, y
  • análisis de puntos extremos

Descriptiva de una variable – frecuencias y cuantiles

Tamaño de la muestra o número de datos

length(mtcars$disp)
## [1] 32
n <- sum(!is.na(mtcars$disp))
n
## [1] 32

Frecuencias absolutas y relativas

Para una variable cuantitativa (numeric)…

frec_abs <- table(cut(mtcars$disp,
      breaks=c(0,100,200,300,400,500,Inf)))
frec_abs
## 
##   (0,100] (100,200] (200,300] (300,400] (400,500] (500,Inf] 
##         5        11         5         8         3         0
frec_rel <- frec_abs / n
frec_rel
## 
##   (0,100] (100,200] (200,300] (300,400] (400,500] (500,Inf] 
##   0.15625   0.34375   0.15625   0.25000   0.09375   0.00000

Cuantiles

quantile(mtcars$disp,.1)
##   10% 
## 80.61
quantile(mtcars$disp,1:4*.2)
##    20%    40%    60%    80% 
## 120.14 160.00 275.80 350.80

Descriptiva de una variable – tendencia central

mean(mtcars$disp)
## [1] 230.7219
median(mtcars$disp)
## [1] 196.3

Descriptiva de una variable – dispersión

var(mtcars$disp) # el denominador es n-1
## [1] 15360.8
sd(mtcars$disp)
## [1] 123.9387
range(mtcars$disp)
## [1]  71.1 472.0
diff(range(mtcars$disp))
## [1] 400.9
IQR(mtcars$disp)
## [1] 205.175
mad(mtcars$disp) # mean absolute deviation
## [1] 140.4764

Descriptiva de una variable – puntos extremos

lmin <- quantile(mtcars$disp,.25) - 1.5*IQR(mtcars$disp) 
lmax <- quantile(mtcars$disp,.75) + 1.5*IQR(mtcars$disp) 
mtcars$disp[mtcars$disp > lmax | mtcars$disp < lmin]
## numeric(0)
lmin <- quantile(mtcars$qsec,.25) - 1.5*IQR(mtcars$qsec) 
lmax <- quantile(mtcars$qsec,.75) + 1.5*IQR(mtcars$qsec) 
mtcars$qsec[mtcars$qsec > lmax | mtcars$qsec < lmin]
## [1] 22.9

Descriptiva de una variable – gráficos

Gráficamente, son útiles para describir un conjunto de datos

  • para variables cuantitativas: el diagrama de caja y el histograma
boxplot(mtcars$qsec)
hist(mtcars$qsec)
  • para variables cualitativas: el diagrama de barras
barplot(table(factor(mtcars$cyl)))

Descriptiva de dos variables

Para el estudio y descripción de la relación entre dos variables, las técnicas más habituales incluyen

  • las tablas de contingencia,
  • los coeficientes de correlación, y muy especialmente el producto-momento de Pearson, y
  • los gráficos de dispersión

Descriptiva de dos variables – tabla de contingencia

Es el concepto equivalente a las tablas de frecuencias en el caso univariante.

Para variables cualitativas…

table(factor(mtcars$am),factor(mtcars$cyl))
##    
##      4  6  8
##   0  3  4 12
##   1  8  3  2

Para variables cuantitativas…

table(cut(mtcars$disp,breaks = seq(0,500,by=100)),
      cut(mtcars$qsec,breaks = seq(13,25,by=3)))
##            
##             (13,16] (16,19] (19,22] (22,25]
##   (0,100]         0       3       2       0
##   (100,200]       1       7       2       1
##   (200,300]       0       3       2       0
##   (300,400]       4       4       0       0
##   (400,500]       0       3       0       0

Descriptiva de dos variables – coeficiente de correlación

Producto-momento de Pearson

cor(mtcars$disp,mtcars$qsec)
## [1] -0.4336979

Coeficiente de correlación de Spearman

cor(mtcars$disp,mtcars$qsec, method = "spearman")
## [1] -0.4597818

Descriptiva de dos variables – gráfico de dispersión

plot(mtcars$disp,mtcars$qsec)

Distribuciones de probabilidad

La distribución de probabilidad de un número aleatorio corresponde a la abstracción teórica de la densidad de los datos de una conjunto. Representa la densidad que se obtendría cunado se generasen infinitos números aleatorios de acuerdos a un mismo procedimiento o experimento.

Un gran número de distribuciones experimentales tienen modelos teóricos con los que se relacionan. Los tres más comunes son la distribución normal (para datos de variables cuantitativas continuas), la distribución uniforme (para cualquier tipo de datos) y la distribución binomial (para el número de veces que se produce un suceso de una determinada probabilidad).

En R, todas las distribuiciones comparten el mismo sistema de funciones

  • r<dist>: genera aleatorios deacuerdo con una distribución de probabilidad
  • d<dist>: devuelve la densidad de probablidad para un valor de variable
  • q<dist>: devuelve el valor de la variable tal que la probablidad
  • p<dist>: devuelve la probablidad que un valor sea inferior o igual a x

Por ejemplo, para la distribución normal estándar…

rnorm(10)
##  [1] -0.5040593 -0.4025292 -1.2126037 -0.7030878 -0.2707858  0.2570524
##  [7]  0.6062162 -0.1483702 -1.3017443  0.4667901
dnorm(0)
## [1] 0.3989423
qnorm(.95)
## [1] 1.644854
pnorm(1.64)
## [1] 0.9494974

Distribuciones de probabilidad - normal

df <- data.frame(x = rnorm(1000, mean = 3, sd = 1))
dfT <-data.frame(x = seq(0,6,length.out=101),
      y = dnorm(seq(0,6,length.out=101),mean=3,sd=1))
ggplot(df, aes(x = x, y=..density..)) +
  geom_histogram(binwidth = .1,fill="#dddddd",col="black") +
  geom_density() +
  geom_line(data=dfT,aes(x=x,y=y),col="red")+
  theme_bw()

pnorm(5,mean = 3,sd = 1)
## [1] 0.9772499
qnorm(.98,mean = 3,sd = 1)
## [1] 5.053749

Distribuciones de probabilidad - uniforme

df <- data.frame(x = runif(1000, min = 10, max = 20))
dfT <-data.frame(x = seq(10,20,length.out=101),
      y = dunif(seq(10,20,length.out=101), min=10, max=20))
ggplot(df, aes(x = x, y=..density..)) +
  geom_histogram(binwidth = .5,fill="#dddddd",col="black") +
  geom_density() +
  geom_line(data=dfT,aes(x=x,y=y),col="red")+
  theme_bw()

punif(12, min=10, max=20)
## [1] 0.2
qunif(.90, min=10, max=20)
## [1] 19

Distribuciones de probabilidad - binomial

df <- data.frame(x = rbinom(100,5,prob=0.5))
df <- df %>% group_by(x) %>%
  summarise(y = n()/nrow(df))
dfT <- data.frame(x = 0:5,
    y = dbinom(0:5,size=5,prob=0.5))

dfJ <- rbind(cbind(teo=FALSE,df),
             cbind(teo=TRUE,dfT))

ggplot(dfJ,aes(x=x,y=y,fill=teo)) +
  geom_bar(position="dodge", stat="identity",
           col="black") +
  scale_fill_manual(values=c("#dddddd","red")) +
  theme_bw() +
  theme(legend.position = "none")

pbinom(2,5,prob=0.5)
## [1] 0.5
qbinom(.5,5,.5)
## [1] 2

Distribuciones de probabilidad – otras

R incorpora otras muchas distribuciones de probabilidad que pueden consultarse en la página correspondiente de la ayuda.

? "Distributions"

Inferencia

En términos estadísticos, la inferencia consiste en la obtención de información sobre la población –el conjunto de los valores existentes– a partir de una muestra representativa de la misma.

Si los datos que tenemos, no constituyen una muestra representativa cualquier información que queramos extrapolar a la población estará sesgada, es decir, desplazada de su valor real.

Los resultados de la inferencia se suelen mostrar mediante una de las dos formas siguientes:

  • el valor de un parámetro para la población, asociado normalmente a un error o un intervalo de confianza,
  • la probabilidad (p) que unos datos, o cualquiera más extremo que ellos, se produzcan siendo cierta una determinada hipótesis, denominada hipótesis nula.

Presentaremos a continuación, las técnicas más comúnmente usadas para

  • analizar el ajuste a una distribución,
  • estudiar la dispersión de la población, y
  • estimar o comparar valores de tendencia central –localización o posición–

Inferencia – distribución

Para contrastar la hipótesis que un conjunto de datos puede representarse mediante un determinada distribución teórica, las pruebas más comunes son la prueba de bondad de ajuste de chi cuadrado (chisq.test) y la prueba de Kolmogorov-Smirnov (ks.test). En ambos casos, los argumentos son las frecuencias absolutas de cada clase para la distribución experimental y las probabilidades esperadas de acuerdo con la distribución teórica o de referencia.

Ejemplo

Se ha tirado cien veces un dado de 10 caras, y se han obtenido cada una de las caras las veces que se muestran en la tabla siguiente.

1 2 3 4 5 6 7 8 9 10
9 6 7 15 11 11 9 13 9 10

Si el dado fuese normal, todas las caras tendrían la misma probabilidad de aparecer. ¿Lo es?

obs <- c(9, 6, 7, 15, 11, 11, 9, 13, 9, 10)  
exp <- rep(.1, 10)
 
chisq.test(x = obs, p = exp)
## 
##  Chi-squared test for given probabilities
## 
## data:  obs
## X-squared = 6.4, df = 9, p-value = 0.6993

Gráficamente, y aunque sin la misma objetividad, la comparación de distribuciones también suele hacerse mediante el gráfico de cuantiles –en R, es la función qqplot.

Si ambos conjuntos de datos corresponden a la misma distribución, sus cuantiles deben quedar alineados junto a la diagonal del gráfico (especialemente en la región central del mismo).

exp <- rep(1:10, obs)
teo <- rep(1:10, 10)
qqplot(x = exp, y = quantile(teo, (1:length(exp))/length(exp)))
qqline(y = exp, distribution = function(x)  quantile(teo,x),
       probs= c(1/length(exp),1))

Para comprobar si un conjunto de datos puede proceder de una población aleatoria de distribución normal, existen pruebas estadísticas más específicas. Entre las muchas pruebas existentes, son de uso común las pruebas de Shapiro-Wilk y Anderson-Darling.

x1 <- rnorm(20, mean = 4, sd = 5)
x2 <- rbeta(20, shape1 = 5, shape2 = .5, ncp = 4)

shapiro.test(x1)
## 
##  Shapiro-Wilk normality test
## 
## data:  x1
## W = 0.97713, p-value = 0.8919
shapiro.test(x2)
## 
##  Shapiro-Wilk normality test
## 
## data:  x2
## W = 0.59299, p-value = 2.493e-06
if(!require("nortest")) {
  install.packages("nortest")
  library("nortest")
}
ad.test(x1)
## 
##  Anderson-Darling normality test
## 
## data:  x1
## A = 0.24361, p-value = 0.7311
ad.test(x2)
## 
##  Anderson-Darling normality test
## 
## data:  x2
## A = 3.4539, p-value = 5.417e-09

Gráficamente se puede usar el gráfico de cuantiles comentado anteriormente.

qqnorm(x1)
qqline(x1)

qqnorm(x2)
qqline(x2)

Inferencia – dispersión

Las inferencias más habituales respecto a la dispersión de una población constituyen los siguientes casos

  • determinar el intervalo de confianza de la varianza o la desviación típica de la distribución,
  • comparar la dispersión de dos conjuntos de datos,
  • comparar la dispersión de tres o más conjuntos de datos –o comprobar la homocedasticidad de distintos conjuntos de datos–.

Cuando los datos estan normalmente distribuidos, el intervalo de confianza de la varianza se calcula a partir de la distribución de chi cuadrado.

x <- rnorm(50, mean = 0, sd = 2)
df <- length(x) - 1
lower <- var(x) * df / qchisq(1 - 0.05/2, df)
upper <- var(x) * df / qchisq(0.05/2, df)
c(lower = lower, var = var(x), upper = upper)
##    lower      var    upper 
## 3.027845 4.339236 6.738175

Para obtener el intervalo de confianza para la desviación estándar, cuando los datos estan normalmente distribuidos, se calcula la desviación estándar a partir de las varainzas calculadas más arriba.

c(lower = sqrt(lower), sd = sd(x), upper = sqrt(upper))
##    lower       sd    upper 
## 1.740070 2.083083 2.595800

Para la comparación de la dispersión de las poblaciones para dos conjuntos de datos, cuando los datos estan normalmente distribuidos, se lleva a cabo con un prueba de F –función var.test en R–.

x <- rnorm(50, mean = 0, sd = 2)
y <- rnorm(30, mean = 1, sd = 1)
var.test(x, y)
## 
##  F test to compare two variances
## 
## data:  x and y
## F = 4.4128, num df = 49, denom df = 29, p-value = 5.781e-05
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  2.217108 8.302372
## sample estimates:
## ratio of variances 
##            4.41283

Si los datos no estan normalmente distribuidos, la comparación de dispersiones puede hacerse mediante el test de Ansari.

x <- rlnorm(50, meanlog = 2, sdlog = 1)
y <- rlnorm(30, meanlog = 2, sdlog = .2)
ansari.test(x, y)
## 
##  Ansari-Bradley test
## 
## data:  x and y
## AB = 823, p-value = 5.921e-05
## alternative hypothesis: true ratio of scales is not equal to 1

Para comparar las dispersiones de más de un conjunto de datos para los cuáles se pueda asumir normalidad, se usa la prueba de Bartlett –bartlett.test en R–.

x1 <- round(rnorm(20, mean = 1, sd = 2),1)
x2 <- round(rnorm(20, mean = 3, sd = 2),1)
x3 <- round(rnorm(20, mean = 5, sd = 2),1)
list(x1,x2,x3)
bartlett.test(list(x1,x2,x3))
## [[1]]
##  [1]  4.6  3.6  1.1  0.1 -2.8  1.3  3.9  1.7 -2.6  4.0  0.1  4.8  1.5  0.7
## [15] -2.9  1.3  3.1 -3.1  1.4 -0.5
## 
## [[2]]
##  [1]  1.3  3.0  2.8  2.8 -2.0  0.5  4.4  4.2  3.6  7.2  4.2  0.2  0.7  2.9
## [15]  6.5  5.2  4.2  3.6  2.0  4.3
## 
## [[3]]
##  [1] 4.7 3.2 7.7 4.3 6.9 6.0 1.3 4.4 4.0 5.6 4.3 0.9 4.0 8.5 2.2 3.3 5.4
## [18] 6.8 2.9 4.0
## 
##  Bartlett test of homogeneity of variances
## 
## data:  list(x1, x2, x3)
## Bartlett's K-squared = 0.97116, df = 2, p-value = 0.6153

Si los datos no cumplen con la hipótesis de normalidad, en lugar de la prueba de Bartlett se usa la prueba de Levene o la de Brown–Forsythe –leveneTest en el paquete car de R–.

Gráficamente y aunque de forma menos objetiva, la comparación de dispersiones puede hacerse usando gráficos de caja, bien directamente, bien centrando los datos.

x1 <- round(rnorm(20, mean = 1, sd = 2),1)
x2 <- round(rnorm(20, mean = 3, sd = 2),1)
x3 <- round(rnorm(20, mean = 5, sd = 2),1)

xs <- data.frame(
  group = factor(c(rep(1, length(x1)),rep(2, length(x2)),
            rep(3, length(x3)))),
  y = c(x1,x2,x3),
  centered = c(x1 - median(x1),x2 - median(x2),x3 - median(x3)))

ggplot(xs, aes(x=group,y=y)) + geom_boxplot() +
  coord_flip() + theme_bw()
ggplot(xs, aes(x=group,y=centered)) + geom_boxplot() +
  coord_flip() + theme_bw()

Inferencia – localización

Las principales inferencias respecto a la localización –o tendencia central– de los datos corresponden a

  • establecer un intervalo de confianza para la tendencia central de la distribución,
  • comparar la tendencia central de un conjunto de datos con un valor preestablecido,
  • comparar la tendencia central de dos conjuntos de datos,
  • comparar la tendencia central de tres o más conjuntos de datos.

Para establecer el intervalo de confianza de la media, cuando los datos están normalmente distribuidos, se usa la distribución t. El intervalo puede calcularse a partir de la función qt o visualizarlo como resultado de la función t.test.

x <- rnorm(10, mean = 2, sd = .5)
x
##  [1] 1.975816 1.157617 2.268955 3.450344 2.088544 1.991210 1.575897
##  [8] 1.727211 2.054207 2.134882
t.test(x)
## 
##  One Sample t-test
## 
## data:  x
## t = 10.915, df = 9, p-value = 1.718e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  1.619169 2.465768
## sample estimates:
## mean of x 
##  2.042468

Si los datos no están normalmente distribuidos, la función wilcox.test ofrece entre sus resultados el intervalo de confianza para la mediana.

x <- rnorm(10, mean = 3, sd = .5) ^ 3
x
##  [1] 12.84760 44.27621 34.72355 32.56462 35.15335 41.52064 42.10572
##  [8] 25.93625 38.23794 61.32878
wilcox.test(x, conf.int = TRUE)
## 
##  Wilcoxon signed rank test
## 
## data:  x
## V = 55, p-value = 0.001953
## alternative hypothesis: true location is not equal to 0
## 95 percent confidence interval:
##  27.47666 44.27621
## sample estimates:
## (pseudo)median 
##       37.08819

Las dos funciones usadas para obtener los intervalos de confianza permiten comparar la tendencia central de un conjunto de datos con un valor preestablecido. Como se ha indicado antes, la prueba de t require normalidad de los datos; la prueba de Wilcoxon, no.

x <- rnorm(10, mean = 2, sd = .5)
x
##  [1] 1.273417 1.719438 2.449046 2.087944 2.312369 2.219420 1.094407
##  [8] 1.654159 2.895084 2.289877
t.test(x, mu = 1.5)
## 
##  One Sample t-test
## 
## data:  x
## t = 2.8388, df = 9, p-value = 0.01944
## alternative hypothesis: true mean is not equal to 1.5
## 95 percent confidence interval:
##  1.601462 2.397570
## sample estimates:
## mean of x 
##  1.999516
wilcox.test(x ^ 3, mu = 8)
## 
##  Wilcoxon signed rank test
## 
## data:  x^3
## V = 33, p-value = 0.625
## alternative hypothesis: true location is not equal to 8

Para comparar la localización de dos conjuntos de datos, deben tenerse en cuenta los siguientes aspectos

  • ¿los datos de ambas muestras están asociados?
  • ¿los datos estan normalmente distribuidos?
  • ¿tienen ambos conjuntos de datos la misma dispersión?

Si los datos están asociados, es decir que existen pares de datos en ambas muestras tales que comparten fuentes de variación, entonces

  • si ambas muestras están normalmente distribuidas, o la diferencia entre conjuntos de datos está normalmente distribuido, se usa la prueba de t sobre la diferencia de los valores entre ambos conjuntos de datos. En R, puede calcularse la diferencia y luego aplicar la función t.test o bien usar directamente la función t.test con la opción paired = TRUE,
  • si los datos no están normalmente distribuidos se usa la prueba de Wilcoxon –wilcox.test sobre la resta o con la opción paired = TRUE en R–.
x <- rnorm(10, mean = 3, sd = 1)
y <- rnorm(10, mean = 3.5, sd = 1)
list(x = x,y = y)
## $x
##  [1] 3.345022 4.567567 2.504564 2.523306 5.007635 2.739927 3.147583
##  [8] 4.118539 2.384038 3.360721
## 
## $y
##  [1] 1.891288 4.138754 2.288343 5.170342 4.556190 4.928725 5.075908
##  [8] 2.936908 3.863685 6.112612
t.test(x - y)
## 
##  One Sample t-test
## 
## data:  x - y
## t = -1.4082, df = 9, p-value = 0.1927
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.8932433  0.4404728
## sample estimates:
##  mean of x 
## -0.7263852
t.test(x, y, paired = TRUE)
## 
##  Paired t-test
## 
## data:  x and y
## t = -1.4082, df = 9, p-value = 0.1927
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.8932433  0.4404728
## sample estimates:
## mean of the differences 
##              -0.7263852
x <- rlnorm(10, meanlog = 3, sdlog = 1)
y <- rlnorm(10, meanlog = 3.5, sdlog = 1)
list(x = x,y = y)
## $x
##  [1]  16.275599 118.552003   3.979148   4.856603  92.009866  10.225033
##  [7]  18.204673 115.666602   9.336026  18.787804
## 
## $y
##  [1]  22.68449  21.23973  17.02835  99.80393  16.41797  43.08411  16.15939
##  [8]  25.74315  85.43956 165.29406
wilcox.test(x - y)
## 
##  Wilcoxon signed rank test
## 
## data:  x - y
## V = 22, p-value = 0.625
## alternative hypothesis: true location is not equal to 0
wilcox.test(x, y, paired = TRUE)
## 
##  Wilcoxon signed rank test
## 
## data:  x and y
## V = 22, p-value = 0.625
## alternative hypothesis: true location shift is not equal to 0

Si los datos no están asociados, entonces

  • si ambas muestras están normalmente distribuidas, o el conjunto de datos –escalados y centrados en cada muestra– está normalmente distribuido, se usa la prueba de t. En R, se usa la función t.test; la opción var.equal permite escoger entre la prueba calculada con la varianza de los datos, y la aproximación de Welch, que se usa cuando la varianzas no son iguales.
  • si los datos no están normalmente distribuidos se usa la prueba U de Mann-Whitney –wilcox.test en R–.
x <- rnorm(10, mean = 3, sd = 1)
y <- rnorm(14, mean = 3.5, sd = 2)
list(x = x,y = y)
## $x
##  [1] 2.756588 4.963940 2.125890 2.265344 2.790045 2.927530 2.223230
##  [8] 2.508153 3.366607 1.868743
## 
## $y
##  [1] -0.03273222  0.37884725  2.53149830  4.33183571  3.04812589
##  [6]  2.78806269  2.21554483  0.85147586  6.08607800  2.70657752
## [11]  3.09383981  2.51752542  2.33729053  4.52201226
t.test(x, y, var.equal = FALSE)
## 
##  Welch Two Sample t-test
## 
## data:  x and y
## t = 0.21253, df = 20.832, p-value = 0.8338
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.9659723  1.1857600
## sample estimates:
## mean of x mean of y 
##  2.779607  2.669713
x <- rlnorm(10, meanlog = 3, sdlog = 1)
y <- rlnorm(14, meanlog = 3.5, sdlog = 2)
list(x = x,y = y)
## $x
##  [1] 16.420154  6.641600 19.814867  4.027625 25.761614 11.242345 16.520460
##  [8] 51.367157 29.315267 19.482693
## 
## $y
##  [1]  442.209234  185.794834  325.379738   31.431309   99.299656
##  [6]  543.129873 1125.963140   66.345844   33.067248  214.198444
## [11]   54.738450    1.569378  957.049102   99.277018
wilcox.test(x, y)
## 
##  Wilcoxon rank sum test
## 
## data:  x and y
## W = 12, p-value = 0.0002743
## alternative hypothesis: true location shift is not equal to 0

Para la comparación de tres o más grupos, las técnicas más habituales son

  • si los datos de los distintos conjuntos de datos están normalmente distribuidos y presentan dispersiones comparables –homocedasticidad–, la técnica recomendada es el análisis de la varianza –aov en R–;
  • si los datos no están normalmente distribuidos pero las distribuciones tienen formas parecidas y se mantiene la homocedaticidad, se usa la prueba de Kruskal-Wallis –en R, krukal.test–.
x <- c(rnorm(10, mean = 3, sd = 1),
       rnorm(8, mean = 3.2, sd = 1),
       rnorm(10, mean = 4, sd = 1))
g <- factor(c(rep(1,10), rep(2,8), rep(3,10)))
test <- aov(x ~ g)
summary(test)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## g            2  5.094  2.5469   2.824 0.0784 .
## Residuals   25 22.546  0.9019                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
x <- c(runif(10, min = 2, max = 4),
       runif(8, min = 1.5, max = 3.5),
       runif(10, min = 3, max = 5))
g <- factor(c(rep(1,10), rep(2,8), rep(3,10)))
kruskal.test(x ~ g)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  x by g
## Kruskal-Wallis chi-squared = 16.925, df = 2, p-value = 0.0002113

Tanto el análisis de la varianza como la prueba de Kruskal-Wallis solo permiten contrastar si algún par de grupos presentan valores de localización distintos, pero no permiten identificar cuál o cuáles.

Para identificar qué pares de grupos presentan diferencias estadísticamente significativos se usan pruebas post-hoc

  • En el caso del análisi de varianza, se puede usar la prueba de Tukey –TukeyHSD en R,
  • Para la prueba de Kruskal-Wallis, se puede usar la prueba de Dunn –dunn.test en el paquete dunn.test de R–

Gráficamente y aunque con menos objetividad, la comparación de la loclización puede hacerse usando gráficos de caja.

x <- c(rnorm(10, mean = 3, sd = 1),
       rnorm(8, mean = 3.2, sd = 1),
       rnorm(10, mean = 4, sd = 1))
g <- factor(c(rep(1,10), rep(2,8), rep(3,10)))

xs <- data.frame(y = x, group = g)
  
ggplot(xs, aes(x=group,y=y)) + geom_boxplot() +
  coord_flip() + theme_bw()

Ajuste de modelos

El ajuste de modelos se expresa de forma habitual en R mediante un objeto formula.

En esta notación, se indica la relación entre variables mediante una expresión de tres términos donde la variable dependiente se indica a la izquierda y las variables dependientes a la derecha.

form1 <- y ~ x              # recta
class(form1)
## [1] "formula"
form1 <- y ~ log(x)         # logaritmo
form1 <- y ~ poly(x,4)      # polinomio de grado 4
form1 <- y ~ x + 0          # recta que pasa por el origen
form1 <- y ~ I(x^.5)        # raíz cuadrada

Ajuste de modelos – regresión lineal

El ajuste de un modelo lineal por mínimos cuadrados ordinarios se hace en R mediante la función lm.

fit1 <- lm(mpg ~ wt, data=mtcars)
summary(fit1)
## 
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5432 -2.3647 -0.1252  1.4096  6.8727 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
## wt           -5.3445     0.5591  -9.559 1.29e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared:  0.7528, Adjusted R-squared:  0.7446 
## F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

Una vez ajustado el modelo, este puede usarse para predecir nuevos datos (o calcular los valores ajustados) mediante la función predict.

df <- data.frame(x=mtcars$wt, y=mtcars$mpg,
  predict(fit1, newdata=mtcars, interval="prediction"))
head(df)
##                       x    y      fit      lwr      upr
## Mazda RX4         2.620 21.0 23.28261 16.92894 29.63628
## Mazda RX4 Wag     2.875 21.0 21.91977 15.59072 28.24882
## Datsun 710        2.320 22.8 24.88595 18.48644 31.28546
## Hornet 4 Drive    3.215 21.4 20.10265 13.78568 26.41962
## Hornet Sportabout 3.440 18.7 18.90014 12.57806 25.22223
## Valiant           3.460 18.1 18.79325 12.47021 25.11630

El modelo se puede visualizar gráficamente a partir de estas predicciones.

df2 <- data.frame(
  wt = seq(min(mtcars$wt), max(mtcars$wt), length.out = 201),
  predict(fit1, 
  newdata = data.frame(wt = seq(min(mtcars$wt), max(mtcars$wt),
                             length.out = 201)), 
  interval="prediction"))

ggplot(mtcars, aes(x=wt, y=mpg)) +
  geom_ribbon(data=df2, mapping = aes(x=wt, ymin=lwr,ymax=upr),
              inherit.aes = FALSE, fill="#dddddd", alpha=0.8)+
  geom_point(shape=3) +
  geom_line(data=df2, mapping = aes(x=wt, y=fit)) +
  theme_bw()

Ajuste de modelos – comprobación

Para comprobar la adecuación del modelo y del método de ajuste usado (OLS) deben analizarse los residuales.

Estos deben ser aleatorios, estar normalmente distribuidos, mostrar homocedasticidad y no mostrar evidencia de puntos especialmente influyentes.

Estas comprobaciones pueden hacerse a partir de los datos de los residuales y las pruebas estadísticas oportunas.

fit1$residuals
##           Mazda RX4       Mazda RX4 Wag          Datsun 710 
##          -2.2826106          -0.9197704          -2.0859521 
##      Hornet 4 Drive   Hornet Sportabout             Valiant 
##           1.2973499          -0.2001440          -0.6932545 
##          Duster 360           Merc 240D            Merc 230 
##          -3.9053627           4.1637381           2.3499593 
##            Merc 280           Merc 280C          Merc 450SE 
##           0.2998560          -1.1001440           0.8668731 
##          Merc 450SL         Merc 450SLC  Cadillac Fleetwood 
##          -0.0502472          -1.8830236           1.1733496 
## Lincoln Continental   Chrysler Imperial            Fiat 128 
##           2.1032876           5.9810744           6.8727113 
##         Honda Civic      Toyota Corolla       Toyota Corona 
##           1.7461954           6.4219792          -2.6110037 
##    Dodge Challenger         AMC Javelin          Camaro Z28 
##          -2.9725862          -3.7268663          -3.4623553 
##    Pontiac Firebird           Fiat X1-9       Porsche 914-2 
##           2.4643670           0.3564263           0.1520430 
##        Lotus Europa      Ford Pantera L        Ferrari Dino 
##           1.2010593          -4.5431513          -2.7809399 
##       Maserati Bora          Volvo 142E 
##          -3.2053627          -1.0274952

Gráficamente, los mismas comprobaciones pueden hacerse a partir de la representación gráfica del modelo.

plot(fit1, which=1:6)

Estas mismas ideas son la base de la aplicación en R de las técnicas de ajuste multilineal (lm), ajuste lineal generalizado (glm) y ajuste no lineal por mínimos cuadrados (nls).

Importación y exportación de datos

Archivos de datos

A menudo es necesario leer y guardar los datos en algún formato tal que permita la importación y exportación de los mismos y su intercambio con otros programas o entornos.

Aunque puede importar y exportar archivos de muy distintos formatos, solo comentaremos como trabajar con

  • archivos de texto, y
  • archivos RData, el formato de almacenamiento de datos de R.

Para los distintos ejemplos usaremos el conjunto de datos flights del paquete nycflights13.

if(!require("nycflights13")) {
  install.packages("nycflights13")
  library("nycflights13")
}
str(flights)
## Classes 'tbl_df', 'tbl' and 'data.frame':    336776 obs. of  19 variables:
##  $ year          : int  2013 2013 2013 2013 2013 2013 2013 2013 2013 2013 ...
##  $ month         : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ day           : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ dep_time      : int  517 533 542 544 554 554 555 557 557 558 ...
##  $ sched_dep_time: int  515 529 540 545 600 558 600 600 600 600 ...
##  $ dep_delay     : num  2 4 2 -1 -6 -4 -5 -3 -3 -2 ...
##  $ arr_time      : int  830 850 923 1004 812 740 913 709 838 753 ...
##  $ sched_arr_time: int  819 830 850 1022 837 728 854 723 846 745 ...
##  $ arr_delay     : num  11 20 33 -18 -25 12 19 -14 -8 8 ...
##  $ carrier       : chr  "UA" "UA" "AA" "B6" ...
##  $ flight        : int  1545 1714 1141 725 461 1696 507 5708 79 301 ...
##  $ tailnum       : chr  "N14228" "N24211" "N619AA" "N804JB" ...
##  $ origin        : chr  "EWR" "LGA" "JFK" "JFK" ...
##  $ dest          : chr  "IAH" "IAH" "MIA" "BQN" ...
##  $ air_time      : num  227 227 160 183 116 150 158 53 140 138 ...
##  $ distance      : num  1400 1416 1089 1576 762 ...
##  $ hour          : num  5 5 5 5 6 5 6 6 6 6 ...
##  $ minute        : num  15 29 40 45 0 58 0 0 0 0 ...
##  $ time_hour     : POSIXct, format: "2013-01-01 05:00:00" "2013-01-01 05:00:00" ...

Partiremos de un subconjunto de flights, para ello empezamos segmentando el data.frame.

fl_ny2ws <- flights[flights$dest %in% c("IAD","BWI"),
                    c("origin","dest","carrier","arr_delay","air_time")]
head(fl_ny2ws)
## # A tibble: 6 x 5
##   origin  dest carrier arr_delay air_time
##    <chr> <chr>   <chr>     <dbl>    <dbl>
## 1    LGA   IAD      EV       -14       53
## 2    LGA   BWI      WN       -19       40
## 3    EWR   IAD      EV        12       52
## 4    JFK   BWI      MQ       851       41
## 5    JFK   IAD      B6        14       52
## 6    LGA   IAD      EV         5       59

Escribir y leer archivos de texto

Para escribir un archivo de texto (delimitado) desde un data.frame se usa la función write.table o cualquiera de sus derivadas. Para leer, la función es read.table.

write.table(fl_ny2ws, file="fl_ny2ws.csv", sep=",", dec=".",
            quote=TRUE, fileEncoding="UTF-8", row.names=FALSE)
fl_ny2ws <- read.table("fl_ny2ws.csv", sep=",", dec=".",
           quote="\"", fileEncoding="UTF-8", header=TRUE)

Si se desea o dispone de un archivo delimitado por tabulaciones, el carácter tabulador se indica como \t.

Escribir y leer archivos de datos de R

Para leer y guardar datos en el formato propio de R, se usan las funciones save y load. Estas permiten almacenar y recuperar cualquier conjunto de variables del entorno de trabajo. Al recuperarlas se recuperan con el mismo nombre con el que se almacenaron.

save(fl_ny2ws, file="fl_ny2ws.rda")
print(load("fl_ny2ws.rda"))   #print muestra el nombre de los objetos

Manipulación avanzada de tablas de datos

Tabla de datos o data frame

Com se ha viso anteriormente, el data frame es el tipo de dato más usado para almacenar tablas de datos.

A menudo, para poder realizar gráficos y/o aplicar distintos procedimientos estadísticos, es necesario manipular la tabla de datos. A esto dedicaremos este apartado.

Entre las operaciones habituales que haremos sobre una tabla de datos, hay

  • renombrar filas o columnas,
  • añadir columnas o filas,
  • segmentar,
  • cambiar el formato de un conjunto de datos,
  • eliminar filas o columnas,
  • crear tablas de datos de resumen, y
  • unir tablas

Terminaremos este bloque introduciendo los funciones básicas del paquete dplyr que facilita la realización de algunas de estas operaciones.

Partimos de una tabla de datos sintética…

df <- data.frame(1:5, letters[1:5], c(rep("a", 3), rep("b", 2)))
df
##   X1.5 letters.1.5. c.rep..a...3...rep..b...2..
## 1    1            a                           a
## 2    2            b                           a
## 3    3            c                           a
## 4    4            d                           b
## 5    5            e                           b

Renombrar filas o columnas

colnames(df) <- c("var1", "var2", "var3") 
rownames(df) <- paste("subject00", 1:5, sep = "")
df
##            var1 var2 var3
## subject001    1    a    a
## subject002    2    b    a
## subject003    3    c    a
## subject004    4    d    b
## subject005    5    e    b

Añadir columnas o filas

df2 <- cbind(df, rnorm(5)) # añadir un vector al data frame
df2$var5 <- 5:1 # assignando valores a una nueva variable
df2
##            var1 var2 var3  rnorm(5) var5
## subject001    1    a    a -1.304938    5
## subject002    2    b    a  1.708467    4
## subject003    3    c    a -1.855023    3
## subject004    4    d    b  1.526245    2
## subject005    5    e    b -1.995378    1
df2 <- rbind(df, list(6, "e", "b"))
df2
##            var1 var2 var3
## subject001    1    a    a
## subject002    2    b    a
## subject003    3    c    a
## subject004    4    d    b
## subject005    5    e    b
## 6             6    e    b

Segmentar

Existen tres formas básicas de seleccionar filas o columnas de una tablas de datos

  • mediante índices numéricos,
  • usando los nombres de filas y columnas, y
  • mediante vectores lógicos

Segmentar – mediante índices

df[1:3,]
##            var1 var2 var3
## subject001    1    a    a
## subject002    2    b    a
## subject003    3    c    a
df[,c(1,3)]
##            var1 var3
## subject001    1    a
## subject002    2    a
## subject003    3    a
## subject004    4    b
## subject005    5    b
df[-3,-2]
##            var1 var3
## subject001    1    a
## subject002    2    a
## subject004    4    b
## subject005    5    b

Segmentar – usando nombres

df[,"var2"]
## [1] a b c d e
## Levels: a b c d e
df$var3
## [1] a a a b b
## Levels: a b
df[,c("var2","var3")]
##            var2 var3
## subject001    a    a
## subject002    b    a
## subject003    c    a
## subject004    d    b
## subject005    e    b

Segmentar – mediante vectores lógicos

df[c(T,T,F,T,F), c(T,F,T)]
##            var1 var3
## subject001    1    a
## subject002    2    a
## subject004    4    b
df[df[,1] == 3 | df[,3] == "b",]
##            var1 var2 var3
## subject003    3    c    a
## subject004    4    d    b
## subject005    5    e    b

Cambiar el formato de un conjunto de datos

Un mismo conjunto de datos puede representar en una tabla de datos de acuerdo con distintas oragnizaciones o formatos.

Se denomina formato ordenado, tidy, aquella organización en la que las filas representan individuos, las columnas, variables de medidas, y las intersecciones los valores de dichas medidas.

Cuando existen más de una variable que corresponde al mismo tipo de medida, entonces los datos pueden organizarse de dos formas

  • usando una columna distinta para cada variable, formato ancho y tidy, o
  • usando una columna para los valores y otra para indicar cuál es la variable representada en esta fila.

El formato más adecuado dependerá de los análisis y visualizaciones que quieran hacerse.

Partiremos de un subconjunto de flights…

if(!require("nycflights13")) {
  install.packages("nycflights13")
  library("nycflights13")
}
fl_ny2ws_W <- flights[flights$dest %in% c("IAD","BWI"),
                    c("origin","dest","carrier","arr_delay","dep_delay")]
fl_ny2ws_W <- cbind(key = 1:nrow(fl_ny2ws_W), fl_ny2ws_W)
head(fl_ny2ws_W)
##   key origin dest carrier arr_delay dep_delay
## 1   1    LGA  IAD      EV       -14        -3
## 2   2    LGA  BWI      WN       -19        -1
## 3   3    EWR  IAD      EV        12        24
## 4   4    JFK  BWI      MQ       851       853
## 5   5    JFK  IAD      B6        14        15
## 6   6    LGA  IAD      EV         5        10

Cambiar el formato de un conjunto de datos – ancho a largo

fl_ny2ws_L <- rbind(
  cbind(edge = rep("origin", nrow(fl_ny2ws_W)), fl_ny2ws_W[,c(1,4)],
        airport = fl_ny2ws_W[,2], delay = fl_ny2ws_W[,6]),
  cbind(edge = rep("dest", nrow(fl_ny2ws_W)), fl_ny2ws_W[,c(1,4)],
        airport = fl_ny2ws_W[,3], delay = fl_ny2ws_W[,5]))
colnames(fl_ny2ws_L) <- c("edge","key","carrier","airport","delay")
head(fl_ny2ws_L)
##     edge key carrier airport delay
## 1 origin   1      EV     LGA    -3
## 2 origin   2      WN     LGA    -1
## 3 origin   3      EV     EWR    24
## 4 origin   4      MQ     JFK   853
## 5 origin   5      B6     JFK    15
## 6 origin   6      EV     LGA    10
tail(fl_ny2ws_L)
##       edge  key carrier airport delay
## 14957 dest 7476      B6     IAD    -4
## 14958 dest 7477      MQ     BWI    -9
## 14959 dest 7478      EV     IAD   -16
## 14960 dest 7479      EV     IAD   -19
## 14961 dest 7480      9E     IAD   -40
## 14962 dest 7481      9E     BWI     2

Cambiar el formato de un conjunto de datos – largo a ancho

fl_ny2ws_W2p1 <- fl_ny2ws_L[fl_ny2ws_L$edge=="origin",]
fl_ny2ws_W2p2 <- fl_ny2ws_L[fl_ny2ws_L$edge=="dest",]
fl_ny2ws_W2p1 <- fl_ny2ws_W2p1[order(fl_ny2ws_W2p1$key),-1]
fl_ny2ws_W2p2 <- fl_ny2ws_W2p2[order(fl_ny2ws_W2p2$key),-c(1,2)]
fl_ny2ws_W2 <- cbind(fl_ny2ws_W2p1,fl_ny2ws_W2p2[,-1])
colnames(fl_ny2ws_W2) <- c("key", "carrier", "origin", "dep_delay",
                           "dest","arr_delay")
head(fl_ny2ws_W2)
##   key carrier origin dep_delay dest arr_delay
## 1   1      EV    LGA        -3  IAD       -14
## 2   2      WN    LGA        -1  BWI       -19
## 3   3      EV    EWR        24  IAD        12
## 4   4      MQ    JFK       853  BWI       851
## 5   5      B6    JFK        15  IAD        14
## 6   6      EV    LGA        10  IAD         5
tail(fl_ny2ws_W2)
##       key carrier origin dep_delay dest arr_delay
## 7476 7476      B6    JFK        12  IAD        -4
## 7477 7477      MQ    JFK        -5  BWI        -9
## 7478 7478      EV    LGA        -5  IAD       -16
## 7479 7479      EV    JFK        -1  IAD       -19
## 7480 7480      9E    JFK        -7  IAD       -40
## 7481 7481      9E    JFK        38  BWI         2

Eliminar filas o columnas

La forma más habitual de eliminar filas o columnas es segmentando la tabla de datos. Sin embargo, una columna también puede eliminarse asignando la misma a NULL.

df <- data.frame(1:5, letters[1:5], c(rep("a", 3), rep("b", 2)))
colnames(df) <- c("var1", "var2", "var3") 
rownames(df) <- paste("subject00", 1:5, sep = "")
df <- df[-2,]
df$var2 <- NULL
df
##            var1 var3
## subject001    1    a
## subject003    3    a
## subject004    4    b
## subject005    5    b

Si lo que se desea es eliminar una variable del entorno de trabajo entonces se usa la función rm.

rm(df)

Creación de resúmenes a partir de datos en formato ancho

sum_fl_ny2ws <- data.frame(edge=c("origin","dest"))

sum_fl_ny2ws$mean_delay <- apply(fl_ny2ws_W[,c("dep_delay","arr_delay")],2,
                                 mean,na.rm=TRUE)
sum_fl_ny2ws$r_delay <- apply(fl_ny2ws_W[,c("dep_delay","arr_delay")],2,
                      function(x) diff(range(x,na.rm=TRUE)))
sum_fl_ny2ws$s_delay <- apply(fl_ny2ws_W[,c("dep_delay","arr_delay")],2,
                              sd,na.rm=TRUE)
sum_fl_ny2ws
##     edge mean_delay r_delay  s_delay
## 1 origin   16.84267     885 47.85548
## 2   dest   13.11556     904 51.52259

Creación de una tabla de resumen a partir de datos en formato largo

sum_fl_ny2ws <- data.frame(edge=c("origin","dest"))
sum_fl_ny2ws$mean_delay <- by(fl_ny2ws_L$delay,fl_ny2ws_L$edge,
                              mean,na.rm=TRUE)
sum_fl_ny2ws$r_delay <- by(fl_ny2ws_L$delay,fl_ny2ws_L$edge,
                           function(x) diff(range(x,na.rm=TRUE)))
sum_fl_ny2ws$s_delay <- by(fl_ny2ws_L$delay,fl_ny2ws_L$edge,sd,na.rm=TRUE)
sum_fl_ny2ws
##     edge mean_delay r_delay  s_delay
## 1 origin   16.84267     885 47.85548
## 2   dest   13.11556     904 51.52259

dplyr

Estas mismas operaciones que se han comentado más arriba, se pueden llevar a cabo también a partir de un sistema coherente de métodos creado como una gramática para la manipulación de datos, el paquete dplyr –que forma parte de tidyverse–.

http://dplyr.tidyverse.org/

En dplyr las operaciones básicas, se realizan mediante cuatro métodos:

  • select: segmentar columnas,
  • filter: segmentar filas,
  • mutate: crear nuevas columnas, y
  • arrange: ordenar.

Estas se encadenan usando el operador pipe, %>%.

df <- flights %>% select(origin, dest, arr_delay) %>% 
  filter(origin == "LGA" & (dest == "IAD" | dest == "BWI")) %>%
  mutate(arr_delay_h=arr_delay/60) %>% 
  arrange(-arr_delay_h)
df
## # A tibble: 1,818 x 4
##    origin  dest arr_delay arr_delay_h
##     <chr> <chr>     <dbl>       <dbl>
##  1    LGA   IAD       398    6.633333
##  2    LGA   IAD       326    5.433333
##  3    LGA   IAD       312    5.200000
##  4    LGA   IAD       306    5.100000
##  5    LGA   IAD       292    4.866667
##  6    LGA   IAD       281    4.683333
##  7    LGA   IAD       280    4.666667
##  8    LGA   IAD       265    4.416667
##  9    LGA   IAD       262    4.366667
## 10    LGA   IAD       253    4.216667
## # ... with 1,808 more rows

Para la creación de resúmenes a partir de tablas en formato largo, es muy útil y cómoda la combinación group_by y summarize.

df %>% group_by(dest) %>% summarise(mean_delay = mean(arr_delay, na.rm=TRUE))
## # A tibble: 2 x 2
##    dest mean_delay
##   <chr>      <dbl>
## 1   BWI  -7.866667
## 2   IAD  13.708258

Por último, las conversiones entre formatos también pueden hacerse a partir de las funciones spread y gather del paquete tidyr –incluido en tidyverse y que se integra de forma natural en la gramática propuestas por dplyr–.